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Abstract 

Tsunamis are often generated by a moving sea bottom. This pa- 
per deals with the case where the tsunami source is an earthquake. 
The hnearized water-wave equations are solved analytically for vari- 
ous sea bottom motions. Numerical results based on the analytical 
solutions are shown for the free-surface profiles, the horizontal and 
vertical velocities as well as the bottom pressure. 
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1 Introduction 

Waves at the surface of a liquid can be generated by various mechanisms: 
wind blowing on the free surface, wavemaker, moving disturbance on the 
bottom or the surface, or even inside the liquid, fall of an object into the 
liquid, liquid inside a moving container, etc. In this paper, we concentrate 
on the case where the waves are created by a given motion of the bottom. 
One example is the generation of tsunamis by a sudden seafioor deformation. 

There are different natural phenomena that can lead to a tsunami. For 
example, one can mention submarine slumps, slides, volcanic explosions, etc. 
In this article we use a submarine faulting generation mechanism as tsunami 
source. The resulting waves have some well-known features. For example, 
characteristic wavelengths are large and wave amplitudes are small compared 
with water depth. 

Two factors are usually necessary for an accurate modelling of tsunamis: 
information on the magnitude and distribution of the displacements caused 
by the earthquake, and a model of surface gravity waves generation resulting 
from this motion of the seafioor. Most studies of tsunami generation assume 
that the initial free-surface deformation is equal to the vertical displacement 
of the ocean bottom. The details of wave motion are neglected during the 
time that the source operates. While this is often justified because the earth- 
quake rupture occurs very rapidly, there are some specific cases where the 
time scale of the bottom deformation may become an important factor. This 
was emphasized for example by Trifunac and Todorovska [1], who consid- 
ered the generation of tsunamis by a slowly spreading uplift of the seafioor 
and were able to explain some observations. During the 26 December 2004 
Sumatra- Andaman event, there was in the northern extent of the source a 
relatively slow faulting motion that led to significant vertical bottom motion 
but left little record in the seismic data. It is interesting to point out that it 
is the inversion of tide-gauge data from Paradip, the northernmost of the In- 
dian east-coast stations, that led Neetu et al. [2] to conclude that the source 
length was greater by roughly 30% than the initial estimate of Lay et al. [3]. 
Incidentally, the generation time is also longer for landslide tsunamis. 

Our study is restricted to the water region where the incompressible Euler 
equations for potential flow can be linearized. The wave propagation away 
from the source can be investigated by shallow water models which may 
or may not take into account nonlinear effects and frequency dispersion. 
Such models include the Korteweg-de Vries equation [4] for unidirectional 
propagation, nonlinear shallow-water equations and Boussinesq-type models 
[5, 6, 7]. 

Several authors have modeled the incompressible fluid layer as a special 
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case of an elastic medium [8, 9, 10, 11, 12]. In our opinion it may be conve- 
nient to model the liquid by an elastic material from a mathematical point 
of view, but it is questionable from a physical point of view. The crust 
was modeled as an elastic isotropic half-space. This assumption will also be 
adopted in the present study. 

The problem of tsunami generation has been considered by a number of 
authors: see for example [13, 14, 15]. The models discussed in these papers 
lack flexibility in terms of modelling the source due to the earthquake. The 
present paper provides some extensions. A good review on the subject is 
[16]. 

Here we essentially follow the framework proposed by Hammack [17] and 
others. The tsunami generation problem is reduced to a Cauchy-Poisson 
boundary value problem in a region of constant depth. The main extensions 
given in the present paper consist in three-dimensional modelling and more 
realistic source models. This approach was followed recently in [1, 18], where 
the mathematical model was the same as in [17] but the source was different. 

Most analytical studies of linearized wave motion use integral transform 
methods. The complexity of the integral solutions forced many authors [9, 
19] to use asymptotic methods such as the method of stationary phase to 
estimate the far-field behaviour of the solutions. In the present study we 
have also obtained asymptotic formulas for integral solutions. They are useful 
from a qualitative point of view, but in practice it is better to use numerical 
integration formulas [20] that take into account the oscillatory nature of the 
integrals. All the numerical results presented in this paper were obtained in 
this manner. 

One should use asymptotic solutions with caution since they approximate 
exact solutions of the linearized problem. The relative importance of linear 
and nonlinear effects can be measured by the Stokes (or Ursell) number [21]: 

a/h a 

where /c is a wave number, a a typical wave amplitude and h the water depth. 
For ^ 1, the nonlinear effects control wave propagation and only nonlinear 
models are applicable. Ursell [21] proved that near the wave front U behaves 
like 

u ~ ti 

Hence, regardless of how small nonlinear effects are initially, they will become 
important. 

Section 2 provides a description of the tsunami source when the source is 
an earthquake. In Section 3, we review the water-wave equations and provide 
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the analytical solution to the linearized problem in the fluid domain. Section 
4 is devoted to numerical results based on the analytical solution. 

2 Source model 

The inversion of seismic wave data allows the reconstruction of perma- 
nent deformations of the sea bottom following earthquakes. In spite of the 
complexity of the seismic source and of the internal structure of the earth, sci- 
entists have been relatively successful in using simple models for the source. 
One of these models is Okada's model [22]. Its description follows. 

The fracture zones, along which the foci of earthquakes are to be found, 
have been described in various papers. For example, it has been suggested 
that Volterra's theory of dislocations might be the proper tool for a quantita- 
tive description of these fracture zones [23]. This suggestion was made for the 
following reason. If the mechanism involved in earthquakes and the fracture 
zones is indeed one of fracture, discontinuities in the displacement compo- 
nents across the fractured surface will exist. As dislocation theory may be 
described as that part of the theory of elasticity dealing with surfaces across 
which the displacement fleld is discontinuous, the suggestion makes sense. 

As is often done in mathematical physics, it is necessary for simplicity's 
sake to make some assumptions. Here we neglect the curvature of the earth, 
its gravity, temperature, magnetism, non-homogeneity, and consider a semi- 
infinite medium, which is homogeneous and isotropic. We further assume 
that the laws of classical linear elasticity theory hold. 

Several studies showed that the effect of earth curvature is negligible for 
shallow events at distances of less than 20° [24, 25, 26]. The sensitivity 
to earth topography, homogeneity, isotropy and half-space assumptions was 
studied and discussed recently [27]. A commercially available code, ABA- 
CUS, which is based on a finite element model (FEM), was used. Six FEMs 
were constructed to test the sensitivity of deformation predictions to each 
assumption. The author came to the conclusion that the vertical layering of 
lateral inhomogeneity can sometimes cause considerable effects on the defor- 
mation fields. 

The usual boundary conditions for dealing with earth problems require 
that the surface of the elastic medium (the earth) shall be free from forces. 
The resulting mixed boundary- value problem was solved a century ago [28]. 
Later, Steketee proposed an alternative method to solve this problem using 
Green's functions [23]. 
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2.1 Volterra's theory of dislocations 

In order to introduce the concept of dislocation and for simplicity's sake, 
this section is devoted to the case of an entire elastic space, as was done in 
the original paper by Volterra [28]. 

Let O be the origin of a Cartesian coordinate system in an infinite elastic 
medium, Xi the Cartesian coordinates {i = 1,2,3), and e,, a unit vector in 
the positive direction. A force F = Fsk at O generates a displacement 
field u^{P, O) at point P, which is determined by the well-known Somigliana 
tensor 



Ui{P,0) = - — {6ikr^nn- ar^k), with a 



i-Kfi ' ' A + 2/i 

In this relation 6ik is the Kronecker delta, A and /x are Lame's constants, 
and r is the distance from P to O. The coefficient a can be rewritten as 
a = 1/2(1 — z/), where u is Poisson's ratio. Later we will also use Young's 
modulus E, which is defined as 

^ /i(3A + 2/i) 
A + /i 

The notation r^, means dr/dxi and the summation convention applies. 

The stresses due to the displacement field (1) are easily computed from 
Hooke's law: 

aij = XSijUk^k + Kuij + Uj.i). (2) 

One finds 

k / r-i ^\ otF I 3x^x jX^ u S^iXj ~\~ Sfi^jXi 3ijX}^ 



An \ A + /i 

The components of the force per unit area on a surface element are denoted 
as follows: 

rpk k , , 

where the u/s are the components of the normal to the surface element. A 
Volterra dislocation is defined as a surface S in the elastic medium across 
which there is a discontinuity Auj in the displacement fields of the type 

Aui = uf - = Ui + VtijXj, (3) 
^jj. (4) 

Equation (3) in which Ui and VLij are constants is the well-known Weingarten 
relation which states that the discontinuity Awj should be of the type of a 
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rigid body displacement, thereby maintaining continuity of the components 
of stress and strain across S. 

The displacement field in an infinite elastic medium due to the dislocation 
is then determined by Volterra's formula [28] 

MQ) = j jj Au^TtdS. (5) 

Once the surface E is given, the dislocation is essentially determined by 
the six constants Ui and Qij. Therefore we also write 

MQ) = j^ JJ 4{P,Q)iyjdS + ^ JJ{x,aUP,Q)-x.a^i{P,Q)}uidS, (6) 

E E 

where Qij takes only the values Q^, ^23, ^31- Following Volterra [28] and 
Love [29] we call each of the six integrals in (6) an elementary dislocation. 

It is clear from (5) and (6) that the computation of the displacement field 
Uk{Q) is performed as follows. A force Fek is applied at Q, and the stresses 
<jfj{P, Q) that this force generates are computed at the points P{xi) on E. In 
particular the components of the force on E are computed. After multipli- 
cation with prescribed weights of magnitude Awj these forces are integrated 
over E to give the displacement component in Q due to the dislocation on 
E. 



2.2 Dislocations in elastic half-space 

When the case of an elastic half-space is considered, equation (5) remains 
valid, but we have to replace af^ in by another tensor uj^y This can 
be explained by the fact that the elementary solutions for a half-space are 
different from Somigliana solution (1). 

The cufj can be obtained from the displacements corresponding to nuclei 
of strain in a half-space through relation (2). Steketee showed a method of 
obtaining the six uo^^ fields by using a Green's function and derived uj\2i which 
is relevant to a vertical strike-slip fault (see below). Maruyama derived the 
remaining five functions [30]. 

It is interesting to mention here that historically these solutions were first 
derived in a straightforward manner by Mindlin [31, 32], who gave explicit 
expressions of the displacement and stress fields for half-space nuclei of strain 
consisting of single forces with and without moment. It is only necessary to 
write the single force results since the other forms can be obtained by taking 
appropriate derivatives. The method consists in finding the displacement 
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field in Westergaard's form of the Galerkin vector [33]. This vector is then 
determined by taking a linear combination of some biharmonic elementary 
solutions. The coefficients are chosen to satisfy boundary and equilibrium 
conditions. These solutions were also derived by Press in a slightly different 
manner [34]. 



Figure 1: Coordinate system adopted in this study and geometry of the 
source model 

Here, we take the Cartesian coordinate system shown in Figure 1. The 
elastic medium occupies the region < and the xi— axis is taken to 
be parallel to the strike direction of the fault. In this coordinate system, 
u^Xi, X2, X3; ,^1, .^25 (,3) is the ith component of the displacement at (xi, X2, X3) 
due to the jth direction point force of magnitude F at {^i,C,2,^3)- It can be 
expressed as follows [22, 31, 34, 35]: 



X3 = —d 



O 




M^(xi,X2,X3) 



X2, -Xs) - m]^(xi, X2, X3) 
+ulj^{Xi,X2, X3) + X-iUlciXi, X2, X3), 



(7) 
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where 
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"3;3 — ^3 and R^ 



In these expressions -Ri = xi 
Rf + R2 + R^- 

The first term in equation (7), X2, — Xs), is the well-known Somigliana 

tensor, which represents the displacement field due to a single force placed 

(^15^25^3) in an infinite medium [29]. The second term also looks like a 
Somigliana tensor. This term corresponds to a contribution from an image 
source of the given point force placed at (^1, .^2, — ^3) in the infinite medium. 
The third term, ul^{xi,X2, X3), and ul(j{xi,X2, X3) in the fourth term are nat- 
urally depth dependent. When 2:3 is set equal to zero in equation (7), the 
first and the second terms cancel each other, and the fourth term vanishes. 
The remaining term, X2, 0), reduces to the formula for the surface 

displacement field due to a point force in a half-space [22]: 
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In these formulas i?^ 
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In order to obtain the displacements due to the dislocation we need to 
calculate the corresponding ^^-derivatives of the point force solution (7) and 
to insert them in Volterra's formula (5) 



Ui 



1 



Uk dS. 



The ^fc-derivatives are expressed as follows: 
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9uIa 
d^k 

dik 



dik 



F 
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2.3 Finite rectangular source 

Let us now consider a more practical problem. We define the elemen- 
tary dislocations Ui, U2 and U3, corresponding to the strike-slip, dip-slip and 
tensile components of an arbitrary dislocation. In Figure 1 each vector rep- 
resents the direction of the elementary faults. The vector D is the so-called 
Burger's vector, which shows how both sides of the fault are spread out: 
D = u+ - u . 
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A general dislocation can be determined by three angles: the dip angle 
6 of the fault (0 < 5 < tt), the slip or rake angle ^ (0 < ^ < tt), and 
the angle between the fault plane and Burger's vector D. When dealing 
with a geophysical application, an additional angle, the azimuth or strike, 
is introduced in order to provide an orientation of the fault. The general 
situation is schematically described in Figure 2. 



O 



Free surface 



X 




Figure 2: Geometry of the source model and orientation of Burger's vector 
D 

For a finite rectangular fault with length L and width W occurring at 
depth d (Figure 2), the deformation field can be evaluated analytically by 
a change of variables and by integrating over the rectangle. This was done 
by several authors [22, 35, 36, 37, 38]. Here we give the results of their 
computations. The final results are represented below in compact form, using 
Chinnery's notation || to represent the substitution 

= f{x,p) - f{x,p- W) - fix -L,p) + fix -L,p-W), 

where p = y cos 6 + d sin 6. Next we introduce the notation 

q = ysm6 — d cos 6, y = r] cos 5 + g sin 5, d = rjsmS — q cos S 

and 

R^ = ^^ + r^^ + q^ = ^^ + y^ + d^ X'=e + q'. 
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The quantities Ui, U2 and U3 are linked to Burger's vector through the 
identities 

f/i = |D| cos0cos6', f/2 = |D| cos0sin6', U3 = \D\sm(j). 
For a strike-slip dislocation, one has 
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For a dip-slip dislocation, one has 
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For a tensile fault dislocation, one has 
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The terms Ji , . . . , /s are given by 
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parameter 


value 


LJl\J cllig,lc U 


1 o 


Fault depth d, km 


25 


Fault length L, km 


220 


Fault width W, km 


90 


t/j, m 


15 


Young modulus E, GPa 


9.5 


Poisson's ratio u 


0.23 



Table 1: Parameter set used in Figures 3, 4, and 5. 



and if cos 5 = 0, 

h -- 
h -- 
h -- 



eg 



2(A + /i) + 



2(A + ^) 
/i 



R + d {R + dy 



\og{R + T]) 



>^ + l^R + d 
ji ^ sin (5 

X + fiR + d 



Figures 3, 4, and 5 show the free-surface deformation due to the three 
elementary dislocations. The values of the parameters are given in Table 1. 



2.4 Curvilinear fault 

In the previous subsection analytical formulas for the free-surface defor- 
mation in the special case of a rectangular fault were given. In fact, Volterra's 
formula (5) allows to evaluate the displacement field that accompanies fault 
events with much more general geometry. The shape of the fault and Burger's 
vector are suggested by seismologists and after numerical integration one can 
obtain the deformation of the seafioor for more general types of events as well. 

Here we will consider the case of a fault whose geometry is described by 
an elliptical arc (see Figure 6). The parametric equations of this surface are 
given by 

c c 

V) = ^, < e < a, y{^, V) = V, --<ti<-, 



z{i,ri) = -{h + d) + -^Jd^. 
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Figure 3: Dimensionless free-surface deformation z/a due to dip-slip fault- 
ing: = 0,^ = 7r/2, D = (0, t/2,0). Here a is |D| (15 m in the present 
application). The horizontal distances x and y are expressed in kilometers. 
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Figure 4: Dimensionless free-surface deformation z/a due to strike-slip fault- 
ing: = 0, 6' = 0, D = (t/i,0,0). Here a is |D| (15 m in the present 
application). The horizontal distances x and y are expressed in kilometers. 
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Figure 5: Dimensionless free-surface deformation z/a due to tensile faulting: 
(f) = 7r/2, D = (0, 0, t/s). Here a is |D|. The horizontal distances x and y are 
expressed in kilometers. 
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z 




Figure 6: Geometry of a fault with elliptical shape. 
Then the unit normal to this surface can be easily calculated: 




We also need to compute the coefficients of the ffist fundamental form in 
order to reduce the surface integral in (5) to a double Riemann integral. 
These coefficients are 

and the surface element dS is 



1 ^a'^ + eiP- 



dS = VEG - F2 d^dT] = --^ ^ ^ didT]. 

Q \/ a? — 

Since in the crust the hydrostatic pressure is very large, it is natural to 
impose the condition that D ■ n = 0. The physical meaning of this condition 
is that both sides of the fault slide and do not detach. This condition is 
obviously satisfied if we take Burger's vector as 



= ^1— ^ ^ =,0, ^ 

i4 



e(&2-a2)' ' v^a4 + - a?) 
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parameter 



value 



Depth event d, km 
Ellipse semiminor axis a, km 
Ellipse semimajor axis b, km 
Fault width c, km 
Young modulus E, GPa 
Poisson's ratio z/ 



20 
17 
6 

15 
9.5 
0.23 



Table 2: Parameter set used in Figure 7. 



It is evident that D = |D|. 

The numerical integration was performed using a 9-point two-dimensional 
Gauss-type integration formula. The result is presented on Figure 7. The 
parameter values are given in Table 2. 

The example considered in this subsection may not be physically relevant. 
However it shows how Okada's solution can be extended. For a more precise 
modeling of the faulting event we need to have more information about the 
earthquake source and its related parameters. 

After having reviewed the description of the source, we now switch to the 
deformation of the ocean surface following a submarine earthquake. The tra- 
ditional approach for hydrodynamic modelers is to use elastic models similar 
to the model we just described with the seismic parameters as input in order 
to evaluate the details of the seafioor deformation. Then this deformation is 
translated to the free surface of the ocean and serves as initial condition of 
the evolution problem described in the next section. 

3 Solution in fluid domain 

The fluid domain is supposed to represent the ocean above the fault area. 
Let us consider the fluid domain Q shown in Figure 8. It is bounded above by 
the free surface of the ocean and below by the rigid ocean floor. The domain 
O is unbounded in the horizontal directions x and y, and can be written as 



Initially the fluid is assumed to be at rest and the sea bottom to be horizontal. 
Thus, at time t = 0, the free surface and the sea bottom are defined by z = 
and z = —h, respectively. For time t > the bottom boundary moves in a 
prescribed manner which is given by 



n = R^x [-h + ax,y,t),r]{x,y,t)]. 



z = -h + C{x,y,t). 
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Figure 8: Definition of tlie fluid domain and coordinate system 



3 Solution in fluid domain 
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The displacement of the sea bottom is assumed to have all the properties 
required to compute its Fourier transform in x, y and its Laplace transform in 
t. The resulting deformation of the free surface z = rj^x, y, t) must be found. 
It is also assumed that the fluid is incompressible and the flow is irrotational. 
The latter implies the existence of a velocity potential y, z, t) which 
completely describes this flow. By definition of 0, the fluid velocity vector 
can be expressed as q = V0. Thus, the continuity equation becomes 

V-q = A0 = O, {x,y,z)en. (8) 

The potential (j){x, y, z, t) must also satisfy the following kinematic boundary 
conditions on the free-surface and the solid boundary, respectively: 

^ dv^d^dri^d^dri^ z = vix,y,t), (9) 
dz dt dx dx dy dy' ' ' ' 

dz dt dx dx dy dy' 

Assuming that viscous effects as well as capillary effects can be neglected, 
the dynamic condition to be satisfied on the free surface reads 

^ + ^|V0p + (7r/ = O, z = 7]{x,y,t). (11) 

As described above, the initial conditions are given by 

r]{x,y,0) = and a^,y,0) = 0. (12) 

The significance of the various terms in the equations is more transparent 
when the equations are written in dimensionless variables. The new inde- 
pendent variables are 

X = Kx, y = ny, J = Kz, t = at, 

where k is a wavenumber and a is a typical frequency. Note that here the 
same unit length is used in the horizontal and vertical directions, as opposed 
to shallow-water theory. 

The new dependent variables are 

a a aa 

where a is a characteristic wave amplitude. A dimensionless water depth is 
also introduced: 

h = nh. 
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In dimensionless form, and after dropping the tildes, equations (8-11) become 

A0 = 0, (x, y, z) G f2. 



d(p df] f dcf) drj dcf) dr]\ 

~^ = + ' z = Kar]{x,y,t), 

oz at \ ox ox oy oy J 

d(f) dC fd(f)dC d(f)dC\ 
oz at \ ox ox oy oy J 

d(f) 1 ^ ,0 gn ^ , , 

— + -Ka\\/(p\ + — ?7 = 0, z = Kar]{x,y,t). 
Ot Z 0"^ 

Finding the solution to this problem is quite a difficult task due to the 

nonlinearities and the a priori unknown free surface. In this study we linearize 

the equations and the boundary conditions by taking the limit as Ka ^ 0. 

In fact, the linearized problem can be found by expanding the unknown 

functions as power series of a small parameter e := na. Collecting the lowest 

order terms in e yields the linear approximation. For the sake of convenience, 

we now switch back to the physical variables. The linearized problem in 

dimensional variables reads 

A(j) = 0, {x,y,z) eR"^ X [-h,0], (13) 
^ = ?^, ^ = 0, (14) 



dz dt 
d(p_dC 
dz dt 



(15) 



-^ + g7] = 0, z = 0. (16) 

Combining equations (14) and (16) yields the single free-surface condition 

d'^(j) dd) , , 

"^+^7-^ = 0, ^ = 0. (17) 



dt"^ dz 

This problem will be solved by using the method of integral transforms. 
We apply the Fourier transform in (x, y): 

W] = m£)= [ fix,y)e-^^''^+'yUxdy, 



rM/] = fXx,y) = I fXkJ)e'^'''+'yUkdi, 
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and the Laplace transform in time t: 

+00 

2[g] = g{s) = J g{t)e-'' dt. 


For the combined Fourier and Laplace transforms, the following notation is 
introduced: 

+00 

dSl[F{x,y,t)] = F{kJ, s) = y e-*(*^^+^^) dxdy j F{x,y,t)e-'' dt. 



R2 



After applying the transforms, equations (13), (15) and (17) become 

+ e)4> = 0, 



dz^ 

^{kj,-h,s) = sak,i,s), (19) 

s^'^{kj,0,s) + g^{kj,0,s) = 0. (20) 
The transformed free-surface elevation can be obtained from (16): 

f}{k,i,s) = -^'^{k,i,0,s). (21) 

A general solution of equation (18) is given by 

0(fc, i, z, s) = A{k, £, s) cosh(m^) + B{k, i, s) sinh(m^), (22) 

where m = \/k^ + P. The functions A{k^^^s) and B{k,i,s) can be easily 
found from the boundary conditions (19) and (20): 

gsC{kJ,s) 



A{kJ,s) 
B{k,i,s) 



cosh(m/i)[s2 + gmtanh{mh)] 



m cosh(m/i) [s^ + gm tanh(m/i)] 
From now on, the notation 

Lj = a/ gm tanh.{mh) (23) 

will be used. The graphs of uj{m), uj'{m) and uj"{m) are shown in Figure 9. 
Substituting the expressions for the functions A, B in (22) yields 

0(A;, 2, s) = w ^L/ 9 ^ — r\ ( cosh(m2;) sinh(m2;) j . (24) 



cosh(m/i)(s^ + u;^) \ gm 
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Figure 9: Plot of the frequency u(m) = gmta3ih{mh) and its derivatives 
du/dm, d'^uj/dm?. The acceleration due to gravity g and the water depth h 
have been set equal to 1. 
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3.1 Free- surface elevation 

From (21), the free-surface elevation becomes 



cosh(m/i)(s2 + io'^) 



Inverting the Laplace and Fourier transforms provides the general integral 
solution 

fi+ioo 

M2 fi—ioo 

One can evaluate the Laplace integral in (25) using the convolution theorem: 
It yields 

t 

i{kx+ly) 



R2 

This general solution contains as a special case the solution for an ax- 
isymmetric problem, which we now describe in detail. Assume that the 
initial solid boundary deformation is axisymmetric: 



Cix,y)=C{r), r = v^oj^ + y2_ 

The Fourier transform '^[({x, y)] = ({k, €) of an axisymmetric function is also 
axisymmetric with respect to transformation parameters, i.e. 

C(A;, I) = C(m), m := + i^. 

In the following calculation, we use the notation ijj = arctan(£/A;). One has 

C{k,i) = [[ ((r)e-^('=^+^^) dxdy = ! dcj) I ((r)e-*"(^™'*+^^^^*Vcir = 



IR2 

2-77 oo OO 



j d(t) j rC(r)e-*^™'=°'^('^-^)(ir = / rC{r)dr / (e-^™™^'^ + e*^'^'=°'^'^)(i(/.. 
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Using an integral representation of Bessel functions [39] finally yields 

oo 


It follows that 

If f 777,g*'"''cos((/>-i/') r _ 

V{r,t) = J^,J d^pj drnjil-usinur)am,t-r)dr 



+ 00 t 

1 f Jo(mr 



I — n — Trdm / (1 — uj smujT)((m,t — T)dT. 
ZTT J cosh[mh) J 



The last equation gives the general integral solution of the problem in the 
case of an axisymmetric seabed deformation. Below we no longer make this 
assumption since Okada's solution does not have this property. 

In the present study we consider seabed deformations with the following 
structure: 

ax,y,t) ■.= ax,y)T{t). (26) 

Mathematically we separate the time dependence from the spatial coordi- 
nates. There are two main reasons for doing this. First of all we want to 
be able to invert analytically the Laplace transform. The second reason is 
more fundamental. In fact, dynamic source models are not easily available. 
Okada's solution, which was described in the previous section, provides the 
static sea-bed deformation (q{x, y) and we will consider different time depen- 
dencies T{t) to model the time evolution of the source. Four scenarios will 
be considered: 

1. Instantaneous: Tj(t) = H(t), where H(t) denotes the Heaviside step 
function, 

2. Exponential: 

Te{t) = I ^_^t ^ ^ q' with a > 0, 



3. Trigonometric: T^{t) = H{t - to) + |[1 - cos{7rt/to)]H{to - t), 

4. Linear: 

'0, t < 0, 
Ti{t) = { t/to, 0<t<to, 
1, t > to. 
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t 



Figure 10: Typical graphs of Te(t) and Tdt). Here we have set a = 6.2, 
to = 0.7. 

The typical graphs of Tc{t) and Tg{t) are shown in Figure 10. Inserting (26) 
into (25) yields 

^ fi+ioo 

V{x, y^t) = Y^J[ ^M^^^ / 4x4^"^^ dkdi. (27) 
(27r)^ J J cosh(m/i) Ztti J + 

Clearly, ri{x,y,t) depends continuously on the source ({x,y). Physically 
it means that small variations of ( (in a reasonable space of functions such 
as L^) yield small variations of t]. Mathematically this problem is said to be 
well-posed, and this property is essential for modelling the physical processes, 
since it means that small modifications of the ground motion (for example, 
the error in measurements) do not induce huge modifications of the wave 
patterns. 

Using the special representation (26) of seabed deformation and pre- 
scribed time-dependencies, one can compute analytically the Laplace integral 
in (27). To perform this integration, we first have to compute the Laplace 
transform of Tj e,c,«(^)- The results are 
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Inserting these formulas into the inverse Laplace integral yields 

fi+ioo 

1 r e^'s^Tds) , 

T~ / — 2~"* = cosa;t, 
Zm J + 

fi—ioo 
fi+ioo 

1 /• e^*s2Te(s) , / _„t ^ • A 
T~ / ^ — = — -cosu't smujt], 

/i—ioo 
fi+ioo 

1 /• e''s^TJs) , 7^ 
-as 



27ri J s2 + u;2 2(72-072) 



11— too 



(cosc<jt — cos7t + — to)[cosco'(t — to) + cos7t]) 



1 /■ e'**s2Ti(s) sincjt-if(t-to)sincu(t-to) 



-ds 



2-Ki J s2 + ci;2 ujtQ 

fi—ioo 

The final integral formulas for the free-surface elevations with different 
time dependencies are as follows: 

r],{x,y,t) = 777^// -———cos ut dkdi, 

[2n) J J cosh[mn) 

K2 

-a2 f f^(kJ)e'^''='+^yWe-''' -cosLjt-^smujt\ „ ,^ 

Ve{x,y,t) = ——- / / — — dkdi, 

(27r)2 cosh(m/i) \ + uj^ J 



Vcix,y,t) 



72 ff ({k, i)e'^''^+^y^ 



(27r)2 J J 2(72 - cu2) cosh(m/i) 

(cosci;t — cos7t + H(t — to)[cosc<j(t — to) + cos7t]) dkdi, 



{2'kY j J cosh(m/i) \ ci^ti 
3.2 Velocity field 



In some applications it is important to know not only the free-surface 
elevation but also the velocity field in the fiuid domain. One of the goals of 
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this work is to provide an initial condition for tsunami propagation codes. 
For the time being, tsunami modelers take initial seabed deformations and 
translate them directly to the free surface in order to obtain the initial 
condition ri{x,y,0). Since a priori there is no information on the flow ve- 
locities, they take a zero velocity field as initial condition for the velocity: 
V(/)(x, y, z, 0) = 0. The present computations show that it is indeed a very 
good approximation if the generation time is short. 

In equation (24), we obtained the Fourier transform of the velocity po- 
tential (^(x, y, z, t): 

4>{k, £, z, s) = r^tw^ J^^^ 9^ fcosh(m^) - — sinh(m^)) . (28) 



cosh.(mh)(s'^ + u"^) \ gm 

Let us evaluate the velocity field at an arbitrary level z = j3h with — 1 < 
(3 < In the linear approximation the value /3 = corresponds to the 
free surface while (3 = —1 corresponds to the bottom. Next we introduce 
some notation. The horizontal velocities are denoted by u. The horizontal 
gradient {d/dx,d/dy) is denoted by Vh- The vertical velocity component is 
simply w. The Fourier transform parameters are denoted k = {k,£). 

Taking the Fourier and Laplace transforms of 



u{x,y,t) = Vh(l>{x,y,z,t)\ 



z=f3h 



yields 



u{k,i,s) = —i(f){k, i, Ph, s)'k 



gsC{k,i)T{s) ( 

' cosh(pm/i) smh(pm/i) k. 



cosh.{mh){s'^ + uj"^) \ gm 

Inverting the Fourier and Laplace transforms gives the general formula for 
the horizontal velocities: 

^ li+ioo 

s _ ig ffkC{k,i)cosh{m(3h)e'^'''=^^y'^ 1 f sT(s)e"* 
^i^,y,t) - ^JJ coshimh) J l^T^'^^'^'" 



fl — tOO 



kC(A;,£)sinh(m/5/i)e*(*^^+^^) 1 [ s^i:{s)e'^ 



An^ J J mcosh(m/i) 27ri J s"^ + u"^ 

R2 fi-ioo 



ds (ik. 



After a few computations, one finds the formulas for the time dependen- 
cies Tj, Tg and T;. For simplicity we only give the velocities along the free 
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surface (/3 = 0): 

. _ ig f f k({k,i)e'^''''+^y^ smut 
4vr^ J J cosh[mn) u 

iga ff kC(fc,^)e^('=^+^3') / « . ^ „ 

Ue{x,y,t) = -r^ — T\ 1 / u\ - cos ujt + - sm ut ) dk, 



ig ffkak,i)e 



i(kx+£y) 



AtQix"^ J J up' cosh(m/;,) 

(1 — COSLOt — Hit — to)[l — COSLoit — to)]) c^k. 



Next we determine the vertical component of the velocity w[x^ z, t\ It 
is easy to obtain the Fourier-Laplace transform w{k^ z, s) by differentiating 
(28): 

-(1.0 \ sgl{KtyY[s) (s^ 

wik, I, z, s) = = — — , , , „ ^7 — coshfm^j — m smhfm^ 

oz cosh{mh){s^ + u!-') \g 

Inverting this transform yields 

^ fi+ioo 

wixvzt) - 1 // cosh(m.)C(fc,£) 1 [ s^ns)e^ 

R2 

9_ ff msinh(mz)C(fc,0 ^^(fc^+^^) 1 /" sT{s)e'' 

471"^ y J cosh{mh) 2'Ki J + uj'^ ' 

R2 fi—icx) 

for — /i < 2; < 0. One can easily obtain the expression of the vertical velocity 
at a given vertical level by substituting z = (3h in the expression for w. 

The easiest way to compute the vertical velocity w along the free surface 
is to use the boundary condition (14). Indeed, the expression for w can be 
simply derived by differentiating the known formula for f]i^e,c,i{x^y,t). Note 
that formally the derivative gives the distributions 6{t) and 6{t — to) under 
the integral sign. It is a consequence of the idealized time behaviour (such as 
the instantaneous scenario) and it is a disadvantage of the Laplace transform 
method. In order to avoid these distributions we can consider the solutions 
only for t > and t 7^ to- From a practical point of view there is no restriction 
since for any e > we can set t = e or t = to + ^- For small values of e 
this will give a very good approximation of the solution behaviour at these 
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"critical" instants of time. Under this assumption we give the distribution- 
free expressions for the vertical velocity along the free surface: 

Wi{x,y,t) = -— ^ // '—r-, — T-7 — ujsinujtdk, 



4vr2 J J cosh(m/i) 



We{x,y,t) = ^11 , rr^2\ — TT— TT e °* + ^ COS a;t- -sin ) dk, 
"^^"'^''^ = -A^-jj 2(7^-^^)coshM) (""""^-^""^^ 

R2 

+H{t — to) [uj sin u;(t — to) + 7 sin •yt]) dk, 

wi{x,y,t) = — — - // — — — [cosujt- H{t- to) cos uj{t- to)] dk. 

4^0^ J J cosh(m/i) 



3.3 Pressure on the bottom 

Since tsunameters have one component that measures the pressure at the 
bottom (bottom pressure recorder or simply BPR [40]), it is interesting to 
provide as well the expression pb{x, y, t) for the pressure at the bottom. The 
pressure p{x, y, z, t) can be obtained from Bernoulli's equation, which was 
written explicitly for the free surface in equation (11), but is valid everywhere 
in the fluid: 

After linearization, equation (29) becomes 

^ + gz + P = 0. (30) 
ot p 

Along the bottom, it reduces to 

^ + ^7^ + + ^ = 0, z = -h. (31) 



The time-derivative of the velocity potential is readily available in Fourier 
space. Inverting the Fourier and Laplace transforms and evaluating the re- 
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suiting expression ai z = —h gives for the four time scenarios, respectively, 
— cos cut ak, 



dt {2nyJJ cosh2(m/i) 



dt (27r)2 JJ a^ + LU^ \ a J (27r)2 

, ^ -T. e + - cos cut + - sincut dk, 

J J m[a^ + oj^) \ \aJ \a/ J 



= "7777^ // TTT — —[smut-H{t-to)smu{t-to)\ dk. 

dt to{2nyjJ u;cosh^(m/i) 

The bottom pressure deviation from the hydrostatic pressure is then given 

by 

Pb{x,y,t) = - -pgC. 

^'^ z=-h 

Plots of the bottom pressure will be given in Section 4. 

3.4 Asymptotic analysis of integral solutions 

In this subsection, we apply the method of stationary phase in order to 
estimate the far-field behaviour of the solutions. There is a lot of literature 
on this topic (see for example [41, 42, 43, 44, 45]). This method is a classical 
method in asymptotic analysis. To our knowledge, the stationary phase 
method was first used by Kelvin [46] in the context of linear water-wave 
theory. 

The motivation to obtain asymptotic formulas for integral solutions was 
mainly due to numerical difficulties to calculate the solutions for large values 
of X and y. From equation (25), it is clear that the integrand is highly os- 
cillatory. In order to be able to resolve these oscillations, several discretiza- 
tion points are needed per period. This becomes extremely expensive as 
r = Y^x^ + — oo. The numerical method used in the present study 
is based on a Filon-type quadrature formula [20] and has been adapted to 
double integrals with exp[i{kx + £y)] oscillations. The idea of this method 
consists in interpolating only the amplitude of the integrand at discretiza- 
tion points by some kind of polynomial or spline and then performing exact 
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integration for the oscillating part of the integrand. This method seems to 
be quite efficient. 

Let us ffist obtain an asymptotic representation for integral solutions of 
the general form 

7]{x, y,t) = — ' ' ^. T(m, t) dkd£, m = VWTP. (32) 

4vr^ J J cosh[mn) 

Comparing with equation (27) shows that T{m, t) is in fact 

T(m,t) = -i. / ^^e-rf.. 
2m J + uj'^ 

fi—ioo 

For example, we showed above that for an instantaneous seabed deformation 
T{m,t) = cosujt, where u"^ = gmtanhmh. For the time being, we do not 
specify the time behaviour T(s). 

In equation (32), we switch to polar coordinates m and ip = arctan(£//c): 

oo 27r ^ 

V{x,y,'t) = rr~r\ T{m,t)m dipdm 

4vr^ J J cosh(m/i) 



oo 2n 

47r^ J cosh[mh) J 



where (r, </)) are the polar coordinates of {x,y). In the last expression, the 
phase function is $ = mrcos{(f — Stationary phase points satisfy the 
condition d^/dip = 0, which yields two phases: ipi = (f and ip2 = + tt- An 
approximation to equation (32) is then obtained by applying the method of 
stationary phase to the integral over tp: 

oo 

Vir, 0, t) ^ -jl= I (C(m, ^)e^(f — ) + C(m, ^ + vr)e^(— f )) dm. 

VSn^r J cosh(m/i) V / 



This expression cannot be simplified if we do not make any further hypotheses 
on the function T{m,t). 

Since we are looking for the far field solution behaviour, the details of 
wave formation are not important. Thus we will assume that the initial 
seabed deformation is instantaneous: 
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Inserting this particular function T(m, t) in equation (32) yields 
where 

^ r mC(m, ^) /• ^,(.,+^,,os(^-^)) ^^^^^ 



dip dm. 



cosh(mh) 



oo ^ 27r 

m({m, ip) ^ ^j(_^t+„„.cos(v'-v)) 
cosh(m/i) 



The stationary phase function in these integrals is 

$(m, Tp) = mr cos((/9 — '(/') ± ut, oj'^ijn) = gmtanhmh. 
The points of stationary phase are then obtained from the conditions 

oip om 

The first equation gives two points, ipi = ip and '?/'2 = + tt, as before. The 
second condition yields 

r duo 

-cos{^-^,,) = ^— (33) 

Since duo /dm decreases from y/gh to as m goes from to oo (see Figure 
9), this equation has a unique solution for m if \r/t\ < \fgh. This unique 
solution will be denoted by m*. 

For |r| > t^/gh, there is no stationary phase. It means physically that the 
wave has not yet reached this region. So we can approximately set /i ~ 
and I2 ~ 0. From the positivity of the function duo /dm one can deduce 
that tpi = ip is a stationary phase point only for the integral l2- Similarly, 
■ip2 = + n is a stationary point only for the integral Ii. 
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Let us obtain an asymptotic formula for the first integral: 



oo 



J cosh[mh) \ V fnr I 
^ ^ 

oo ^ 

r y cosh(m/i) 





r I y |u;"(m*)|t cosh.{m*h) 



2n CO;?A^j^gi(a;(m*)t-rnV) 

t V cosh{m*h) 

In this estimate we have used equation (33) evaluated at the stationary phase 
point (m*, '02): 

duj , . 

r = t— . 34 

dm 

Similarly one can obtain an estimate for the integral l2'- 

t V — cosh(m*/i) 

Asymptotic values have been obtained for the integrals. As is easily observed 
from the expressions for Ii and I2, the wave train decays as or 1/r, which 
is equivalent since r and t are connected by relation (34). 



4 Numerical results 

A lot of numerical computations based on the analytical formulas ob- 
tained in the previous sections have been performed. Because of the lack of 
information about the real dynamical characteristics of tsunami sources, we 
cannot really conclude which time dependence gives the best description of 
tsunami generation. At this stage it is still very difficult or even impossible. 

Numerical experiments showed that the largest wave amplitudes with the 
time dependence Tc{t) were obtained for relatively small values of the char- 
acteristic time to- The exponential dependence has shown higher amplitudes 
for relatively longer characteristic times. The instantaneous scenario Tj gives 
at the free surface the initial seabed deformation with a slightly lower am- 
plitude (the factor that we obtained was typically about 0.8 ~ 0.94). The 
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Parameter 


Value 


Young modulus, E, GPa 


y.o 


Poisson ratio, z/ 


U.Z 1 


Fault depth, d, km 


on 


Dip angle, 6, ° 


1 Q 
id 


Strike angle, 6, ° 


nn 

yu 


Normal angle, 0, ° 


U 




fin 


Fault width, W, km 


40 


Burger's vector length, D , m 


15 


Water depth, h, km 


4 


Acceleration due to gravity, g, m/s^ 


9.8 


Wave number, k, l/m 


10-^ 


Angular frequency, u, Hz 


10-2 



Table 3: Physical parameters used in the numerical computations 



water has a high-pass filter effect on the initial solid boundary deformation. 
The linear time dependence Ti{t) showed a linear growth of wave amplitude 
from to also ~ 0.9Co, where Co = max \C{x, y)\. 

(x,3/)gM2 

In this section we provide several plots (Figure 11) of the free-surface 
deformation. For illustration purposes, we have chosen the instantaneous 
seabed deformation since it is the most widely used. The values of the pa- 
rameters used in the computations are given in Table 3. We also give plots of 
the velocity components on the free surface a few seconds (physical) after the 
instantaneous deformation (Figure 12). Finally, plots of the bottom dynamic 
pressure are given in Figure 13. 

From Figure 12 it is clear that the velocity field is really negligible in 
the beginning of wave formation. Numerical computations showed that this 
situation does not change if one takes other time-dependences. 

The main focus of the present paper is the generation of waves by a moving 
bottom. The asymptotic behaviour of various sets of initial data propagating 
in a fluid of uniform depth has been studied in detail by Hammack and Segur 
[47, 48] . In particular, they showed that the behaviours for an initial elevation 
wave and for an initial depression wave are different. 
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Figure 11: Free-surface elevation at t = 0.01,0.6,3,5 in dimensionless time. 
In physical time it corresponds to one second, one minute, five minutes and 
eight minutes and a half after the initial seabed deformation. 
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Figure 12: Components u, v and w of the velocity field computed along the 
free surface at t = 0.01, that is one second after the initial seabed deforma- 
tion. 
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Figure 13: Bottom pressure at t = 0.01,0.6,3,5 in dimensionless time. In 
physical time it corresponds to one second, one minute, five minutes and 
eight minutes and a half after the initial seabed deformation. 
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